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Abstract. Finite density quantum field theories have evaded first principle Monte-Carlo 
simulations due to the notorious sign-problem. The partition function of such theories appears as 
the Fourier transform of the generalised density-of-states, which is the probability distribution of 
the imaginary part of the action. With the advent of Wang-Landau type simulation techniques 
and recent advances [1], the density-of-states can be calculated over many hundreds of orders 
of magnitude. Current research addresses the question whether the achieved precision is high 
enough to reliably extract the finite density partition function, which is exponentially suppressed 
with the volume. In my talk, I review the state-of-play for the high precision calculations of 
the density-of-states as well as the recent progress for obtaining reliable results from highly 
oscillating integrals. I will review recent progress for the Z$ quantum field theory for which 
results can be obtained from the simulation of the dual theory, which appears to free of a sign 
problem. 


1. Introduction 

The properties of strongly interacting matter has sparked many important investigations using 
accelerator experiments and large scale theoretical studies. In an astrophysical context, the 
theory of interactions, QCD, dominates the area around 10 -6 seconds after the Big Bang when 
quark matter confines to colour neutral hadrons. Since QCD captures the impact of the strong 
nuclear forces, it is widely believed that QCD plays an essential role to understand compact 
star matter as it may exist nowadays e.g. in neutron stars (see figure [I] for an illustration). The 
so-called QCD phase diagram characterises the states of matter as a function of the density and 
temperature in a 2d graph. Completing this diagram in its extreme regions and in particular 
for cold and dense matter is an outstanding problem which still triggers model building and 
new techniques for computer simulation after 30 years of intense research. Under moderate 
conditions, quarks and gluons are confined to colour-neutral hadrons, while this confinement 
feature of QCD is believed to cease to exist under extreme conditions, temperatures and/or 
densities. In the early universe before the QCD phase transition roughly at time 10 -6 seconds, 
matter has not yet clustered due to gravity and, as a result, the density is very low compared to 
e.g. nuclear matter density. Similarly, the conditions in heavy ion collision experiments carried 
out at RHIC (located at Brookhaven National Laboratory (BNL) in Upton, New York) or LHC 
(built by the European Organization for Nuclear Research (CERN)) produce very hot matter 
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Figure 1. QCD impact on the evolution of the early Universe. 
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at low densities. From an experimental point of view, very little is known even at moderate 
densities let alone the cold and dense regime of compact star matters (see figure [2]). 

Lattice gauge simulations offer first principle non-perturbative results with good control over 
the errors. Markov chain Monte-Carlo simulations can be very successfully applied and reliable 
results for the states of matter can be obtained at all temperatures and small baryon chemical 
potentials (shaded area in figure [3j see e.g. [2] for a review). Over the recent years, many new 
models have been developed to describe QCD matter at high densities: quarkyonic matter is 
motivated by the so-called ’t Hooft limit of the hypothetical theory with a large number of 
colours and sees the quarks deconhned inside the Fermi sphere while only the baryons at the 
surface of the Fermi sphere undergo colour confinement—13]. The chiral magnetic effect [4] 
takes into account the strong magnetic field that are produced by the colliding charges in heavy 
ion experiments and uses an effective quark model to estimate their impact on the QCD phase 
structure. In the confinement phase, quarks in a certain gluonic background can change their 
statistics from Fermi to Bose type. At moderate densities, these so-called centre dressed and 
confined quarks therefore might undergo Bose condensation leading to to new state of cold and 
















confined matter [Sj. None of these ideas have yet been tested by first principle calculations. 
At a finite baryon chemical potential, the QCD action acquires an imaginary part, and, hence, 
standard Monte-Carlo techniques cannot be applied for the lattice simulation of QCD at finite 
densities. This has become known as the notorious sign problem. Early attempts pursued 
Monte-Carlo simulations with the modulus of the quark determinant and included the phase 
factor to the observable to be calculated. It turns out that the expectation value of the phase 
factor is very small (see |6j for an illustration) showing that the such generated Monte-Carlo 
configurations contain very little information on finite density QCD. The sign problem has 
turned into an overlap problem. 

The recent past has seen a renaissance of ideas targeting dense and cold fermionic matter. 
Some of the methods were proposed some time ago, but have now reached an unprecedented 
level of sophistication. The reweighting approach proposed by Fodor and Katz [7] uses a multi¬ 
parameter action to optimise the overlap. This method is particularly suited for intermediate 
temperatures and might possess a reach that covers the critical endpoint of the QCD phase 
diagram |8j. Langevin simulations of lattice gauge theories avoid the positivity constraint of 
the Gibbs factor, which lies at the heart of Monte-Carlo simulations, and might therefore be 
suitable for finite density simulations |9, 10]. This technique regained a lot of interest when Aarts 
showed that stochastic quantisation can evade the sign problem at least for the relativistic Bose 
gas SHE!. Although the conceptional question whether the approach converges to the correct 
answer m is still under active investigations, the approach is one of the few methods that 
are currently applied to finite density QCD P3j. It might appear that integrating the gluonic 
degrees of freedom before the fermion fields alleviates the sign problem. This could be done 
e.g. in the strong coupling limit m leading to a description of Nuclear Physics suiting lattice 
simulations |16j . It was observed in Id QCD that even integrating part of the gluonic degrees 
of freedom leads to substantial improvements P33- Finally, a reformulation of the theory might 
improve on the sign the problem or remove it altogether. Indeed, theories the dual of which 
are real are then accessible by standard Monte-Carlo techniques. Example for complex action 
spin models that are real upon dualisation are spin model [T5j or the 0(2) model at finite 
densities na sol- These models can be efficiently simulated using worm type algorithms |21j . 
Alternative lattice discretisations [22] and spin blocking techniques in combination with the 
(tensor) renormalisation group approach [23] might be equally successful to eliminate the sign 
problem from Yang-Mills theories. Also not hinging on a real dualisation of the theory is 
the fermion bag approach approach by Chandrasekharan [23] for which the sign problem is 
relegated to finite size fermion bags. This approach has been seen to be very efficient for 
fermion theories with four-fermion coupling such as the Thirring model with massless fermions 
on large lattices [ 25] . 

An efficient alternative to conventional Monte-Carlo simulations is based upon the numerical 
computation of the density of states using the nmlti-canonical algorithm m or a Wang-Landau 
type strategy [26] . A modified version of the Wang-Landau method is the LLR algorithm ]T], 
which is effective for theories with continuous degrees of freedom as opposed to spin models (see 
also [28]). The latter algorithm has been extended from the calculation of the action distribution 
to accessing probability distributions of other extensive quantities such as the SU(2) Polyakov 
line [29]. Furthermore, it has been proposed in [30] to use LLR techniques for a high precision 
calculation of the distribution of the imaginary part of the action. Once this quantity has been 
determined, the partition function of the complex theory can be computed semi-analytically by 
carrying out the Fourier transform of the corresponding probability distribution [30]. Below, we 
will summarise the state of affairs concerning density-of-states methods and the LLR algorithm 
in particular to simulate theories with a sign problem. An overview on selected new methods 
solving the sign problem can be found in the recent review by Aarts m • 


2. Density-of-states approach to complex action systems 

2.1. Density-of-states and the overlap problem, 

Let us consider the partition function Z of a theory of one degree of freedom cj) with a complex 
action: 

Z = J V<f> exp {/3S R [(j)\ + i/xS7[0]} (1) 

where /i is the “chemical potential”, and Sr/j are real and imaginary parts of the action. We 
introduce the generalised density of states 113] by 

Pp{s) = j T>(j> S(s - exp{/3Sii[(/>]} . (2) 


For (3 = 0, Po{s) just counts the number of states with the constraint that the imaginary part 
of the action is given by s. At finite (3, the number count is weighted by the “Gibbs” factor 
exp{0Sn [0]}. Once Pp(s) is known, the task to calculate the partition function boils down to 
evaluate the integral 

Z(f3,n) = J ds Pp(s) exp {ijis}. (3) 

Since the integrand in ([2]) is perfectly real, the difficulties with the sign problem are relegated 
to the 1-dimensional integral Q. In fact, Pp(s) could be estimated by performing a standard 
Monte-Carlo simulation with the Gibbs factor exp{/3 Sr [</>]} and to bin the values for Sj in a 
histogram. The LLR algorithm will provide us with 0 = 0, Pq(s) over hundreds of orders of 
magnitude (see below for an example). The aim of this paper is to discuss the options for a 
calculation of the highly oscillating integral (Hi- 

Let us scope the amount of difficulty that resides with this task. We firstly note that the 
action S) is an extensive quantity S/(^») = Vsj with V the number of degrees of freedom (volume) 
and with the action density S{ of order one. A good qualitative choice (see e.g. m) is given by 


Pp(s) = exp 



2 

leading to Z = J ds e~ s ^ exp{i/rs} oc exp{—^-G} 


For a chemical potential fi of order one, Z is exponentially suppressed with the number of degrees 
of freedom V. On the other hand, Pp(s) is only known numerically and of order one at least for 
small s. For a successful evaluation of the integral in Z ©> any numerical method for obtaining 
Pfs(s) must have the properties 

• exponential error suppression for extensive quantities 

• for the whole domain of support for Sj. 

The LLR algorithm proposed in [I] just delivers that. 


2.2. The Zj, spin model as showcase 

The key question is whether the quality of the result for Pp(s) obtained by the LLR algorithm 
is good enough to admit a reliable calculation of the partition function Z via the integral ([3]). 
The answer to this question is model dependent. The 3-dimensional Z% spin model on a cubic 
lattice at finite density maps onto a real action system upon dualisation and is thus open to 
standard Monte-Carlo simulations. It serves as a first benchmark test in the feasibility study 
for our approach to theories with a sign problem. Here, we will review our findings for the 
3-dimensional case. Degrees of freedom are the centre elements z(x ) that take values from the 
group Z 3 


z(x) <E {1 ,z+,z-} 


z± = (l + iVs)/2. 




Figure 4. Histogram count for N + — N. 
Sj (linear scale); 24 3 lattice, t = 0.17, n 
0.05. 


Figure 5. Same histogram on 
a logarithmic scale with the LLR 
result now reaching beyond AC — 
AC ~ 5, 500. 


The action is given by 

S[z] = rY^\ z xZ* x+v + cc} + + ^x] > V = K ex P(^) ) V = « exp(-/r) . (4) 

Jr,*' a; 

This model is inspired by finite density QCD in the heavy quark limit, and the parameter r is 
reminiscent of the temperature and n reflects the quark hopping parameter m- Apparently, 
the action becomes complex as soon as [i ^ 0. If for a given configuration z{x) the quantity N± 
represents the number of spins on the lattice with z = z±, the imaginary part of the action can 
be written as: 

Sj = —— 7 ^ — z*] = \/3«sinh(//) [N + — A^_] . (5) 

2 * X 

We have performed a standard Monte-Carlo simulation using a 24 3 lattice, k = 0.17 and k = 0.05 
to obtain a histogram for A+ — AC (see |30] for details). The result is shown in figure [4] on a 
linear scale (see figure [5] with the y-axis on a logarithmic scale). Note that we have very little 
“events” with AT+ — AC > 1000. For the calculation of the Fourier transform to obtain the 
partition function Z in ([3]), the tails with |A^+ — AC| 1000 significantly contribute for ^ ~ 1. 
We recover the overlap problem in the light of the density-of-states approach. Our result for 
P(N + — AC) using the LLR method is also shown in the figures [4] and [5j We find a reassuring 
agreement with the standard simulation result and moreover we have obtained the distribution 
P(N+ — AC) for values as large as N + — AC = 5000 and over sixty orders of magnitude. 

3. The partition function from highly oscillating integrals 

3.1. Polynomial fit 

Our task is now to carry out the Fourier transform of the generalised density of states P(s) 
in order to obtain the partition function Z (see ([3])). One advantage of our approach is that 
we can use sophisticated integration techniques, which converge like l/n p , p > 1, where n is 
the number of evaluations of the integrand. Note that any Monte-Carlo integration that sub¬ 
sums the Fourier transform necessarily converges at best like 1 /y/n. Note, however, that even 
the sophisticated integrations techniques would fail for sizeable values of the chemical potential 
without further knowledge of the function P(s). We have so far studied the density-of-states for 





































the theories U( 1), SU( 2 ), SU( 3) and Z 3 and a common feature has been that log P{s) is indeed 
very smooth and, as expected, monotonic functions of its variable s. In the present case, we also 
have the symmetry under reflection P(s) = P(—s). The “smoothness” of P{s) is summarised 
by the fact that the Taylor expansion 

9 

In P(s) = > q = 2,4, 6 , 8 ,... ( 6 ) 

i>0,even 


can produce results that are indistinguishable from the numerical findings for P(s) within error 
bars. Depending on the region of the parameter space /3 = (r,«), q as small as 4 might be 
sufficient. As soon as an acceptable representation of the numerical data in terms of the ansatz 
([ 6 ]) is found, the calculation of the partition function can be performed in a “semi-analytic” way 
using advanced numerical integration techniques: 


Z(g) 


2 


ds P(s ) cos(g s ) 


2 


ds exp 



cos (ns) . 


(7) 


3.2. Asymptotic referencing 

Similar to the scenario in the previous subsection, it might be useful to describe the gross features 
of P(s) by an analytic function, say P asy (s) especially for large values s. Decomposing 

P(s) = P(s)P asy (s), (8) 

the function P(s) might have a moderate range of values although P(s) spans many orders 
of magnitude. For a technical side-remark, we point out that the LLR algorithm [1J can be 
easily adapted to directly produce P(s) for a given choice for P asy (s). The partition function is 
obtained by Fourier transformation: 

Z(g) = FT[P](fi) = j ds P{s ) e * 8 = J dx FT[P](x) FT[P asy ]( M - x) . (9) 

The idea is that FT[P asy ] is analytically available and has already incorporated a good deal of 
the cancellations. For moderate values of k and r, a sensible choice is 

Pasy(s) oc exp{-as 2 } leading to FT[P asy ](/it - x) oc exp |- ^ ^ | . (10) 

Depending on the size of the intrinsic scale a, we only need to numerically calculate P(s) for 
s ~ g for the chemical potential g of interest. 

3.3. Eigenfunction expansion 

Let us expand the generalised density-of-states P(s) in terms of a complete set of eigenfunctions 
(in the L2 sense) that are eigenfunctions of a differential operator: 

P(s ) = ^CnMks), (11) 

n 

where A: is a parameter that will be adapted to the intrinsic scale of the theory under studies. 
We need not necessarily adopt an eigensystem with an all discrete set of eigenfunctions. Here, 
we have indeed the eigenfunctions of the harmonic oscillator in mind having this property: 


-j^if n (ks) + k A s 2 if n (ks) = k 2 (2n + 1) ip n (ks) . 


(12) 



The ortho-normal eigenfunctions are the well-known Hermite functions: 


ipn(x) = 


2 Tl n!y / 7r 


- x2 ' 2 H n (x), H n (x) = {- l) n e* 


dx r 


(13) 


where the H n {x) are the Hermite polynomials. Using the “Schrodinger equation” (12), one 
proves that the eigenfunctions are fixed points of the Fourier transformation: 


FT[^]( M ) = / dsMks)e^ s = 


We therefore find for the partition function 


V2^ .. 




i n ^ m • 


(14) 


Z(fj) = FT[P]( M ) = ]Tc n ^ • 


(15) 


If the chemical potential /r is larger than the intrinsic scale k, we observe that we attain 
small values for Z already from the asymptotic behaviour of the Hermite functions, i.e., 
ipn(x) ~ exp{— x 2 /2}. 


4. The Z 3 spin model - numerical results 



Figure 6 . The overlap 0{/j) from a direct Monte-Carlo simulation of the dual theory 
(DS) and from the LLR approach for several lattice sizes L 3 . t = 0.1, k = 0.01. 


In order to quantify the influence of the imaginary part, we introduce the partition function 
Zmod that features the Z 3 action Q from which we have dropped the imaginary part: 

SmodM = r E[^- z E' + cc \ + K cosh(^) [2JVi + N + + N-] , 

x,u 


(16) 












where N± is the number of z = 1 elements on the lattice. The partition function of the modified 
theory does depend on the chemical potential, but note that it can be simulated by standard 
Monte-Carlo techniques since its Gibbs factor is real and positive by construction. We then 
define the overlap between the full theory and the modified theory by 


O(p) 


Zmod ( m ) 


(17) 


We point out that being able to calculate the overlap O(p) provides access to the observables 
of the full theory. For instance, the density p{p) acquires two parts: 


p(m) 


d In Z(n) 

dp 


d In O(p) d In Z mod (p) 
dp dp 


(18) 


where the latter part is free of a sign problem and is calculable by standard Monte-Carlo 
simulations. 

An appealing feature of the Z 3 spin model is that its dual is a real theory open for Monte-Carlo 
simulations. In fact, the theory can be efficiently simulated by a worm type algorithm (see [18]) 
where the “worms” are conserved flux lines of the dual theory (see ( 20 j for this interpretation). 
It is still difficult to calculate the partition function itself since theories with different chemical 
potentials differ in the free energy density f(p) leading to poor overlap: 


Z(p + A p) 
Z{p) 


exp{- [f{p + A p) - f(p)} VJ . 


(19) 


We used a variant of the “snake algorithm” [32] to calculate ratios of the partition function and 
to reconstruct the partition function from those: 


Z(k A p) 


z{ o) n 


£=1 


Z(£Ap) 
Z((£-l)Ap) ’ 


( 20 ) 


where A p must be chosen small enough (depending on the number of degrees of freedom V ) to 
ensure a good enough signal-to-noise ratio. We here point out an advantage of the LLR approach: 
k simulations of the dual theories are necessary to arrive at the target value pt = kAp, while 
the LLR approach can aim directly at the chemical potential pt of interest. 

We have simulated the Z 3 theory with non-zero chemical potentials using the LLR 
method m- We point out that the method is exact implying that the numerical results need to 
agree within error bars with the exact values. Note, however, that sophisticated techniques for 
the error analysis (e.g. bootstrap) might be needed and that standard Gaussian error analysis 
might fail at least in certain regions of parameter space [33]. Our results from the direct 
simulation (DS) of the dual theory are shown in figure [ 6 ] in comparison with our results from 
the LLR approach. In the latter case, we have used the method of the polynomial fit from 
subsection |3.1| for the evaluation of the highly oscillating integral. We indeed encounter a strong 
sign problem since the overlap is as small as 10 -16 for p ~ 2 . So far, we have only explored a 
limited region of the parameter space (r, k), but our results serve as a proof of concept: for a 
range of volumes and for the case of a strong sign problem, the simulation of the theory in its 
original degrees of freedom is feasible using the LLR techniques. 








5. Conclusions 

Quantum field theory at finite density or, more general, statistical systems with complex 
actions such as the imbalanced Fermi gas [2H still await first principles results from computer 
simulations. In the latter case, theoretical findings can be scrutinised against experiments, 
and, given the level of abstraction that went into model building, agreement would signal an 
understanding of the materials at hand (see e.g. [35]). In the context of QCD at finite baryon 
densities, a variety of mechanisms have been proposed over the last couple of years that should 
describe the states of baryon matter in the intermediate temperature and density range with or 
without a strong magnetic field. Proposals feature “quarkyonic matter” [3], suggested on the 
basis of the ’t Hooft limit, the “chiral magnetic effect” [4j or the “Fermi-Einstein condensation” 
of quarks [5], and this is not meant to be a complete list. Lattice simulations would provide a 
genuine non-perturbative approach with good systematic control of the errors, but are hampered 
by the notorious sign problem: for a non-vanishing chemical potential, the Gibbs factor is 
complex (or, at least, not positive definite) and the action based importance sampling, which is 
at the heart of the Monte-Carlo simulations, are impossible. 

Alongside the new theory proposals for the potential state of matter a finite densities, the last 
decade has seen promising progress for the simulation of theories with a sign problem. In fact, 
many of the related ideas are rooted in the literature for decades, but techniques have reached an 
unprecedented level of sophistication. A good example are the Langevin simulation of complex 
actions systems, which date back to the early works by Parisi [9j and Karsch and Wyld [TOj from 
the mid eighties, but underwent a Renaissance when it was realised the Silver-Blaze problem 
can be avoided for the case of a relativistic Bose gas m- 

Similarly, the LLR algorithm [I] emerged from a modification of Wang-Landau type 
algorithms [233] and have progressed along the lines of the so-called density-of-states methods 
(see e.g. m or [371 ESI sang). The LLR algorithm copes with continuous degrees of freedom 
and is designed to numerically calculate the probability distribution of an extensive quantity, 
the action [lj or e.g. the Polyakov line [29] . to an unprecedented precision and for regions of 
the variable (e.g. action) that would never be visited by an action based importance sampling 
Monte-Carlo approach. Thus, the LLR algorithm naturally solves overlap problems. Given the 
belief of a correspondence between the overlap and the sign problem in finite density quantum 
field theory, it was natural to explore its readiness for complex action theories m- Here, the 
LLR algorithm provides a high quality probability distribution for the imaginary part of the 
action, and the partition function emerges as the Fourier transform of this distribution with 
the chemical potential as its frequency (see (|3]) ). Details and advances of the LLR methods 
have e.g. been reported in flj '28., 33J [30]. In this paper, we have focused on possible techniques 
to extract an signal, which is exponentially small with the volume, from the highly oscillating 
Fourier integral. As a proof of concept, we have studied the Z 3 spin model m- For this model, 
we have used (for a limited range of the parameter space) the Polynomial Fit technique from 
subsection 3.1. Despite of a severe sign problem (as quantified by a phase factor expectation 
value at the 10 16 level; see figure [b]), we were able to obtain reliable results by simulating the 
theory in its original formulation using the LLR techniques. An analysis of the full parameter 
space of the Z 3 model, the LLR simulation of more involved theories (e.g. the 0(2)-model) and 
the exploration of the techniques outlined in section [3] to carry out the Fourier transform are 
currently work in progress. 
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